if (simGetScriptExecutionCount()==0) then
    -- Global variables
    population_size = (20)                  -- the initial population size  default: 20
    max_evaluations = (200)                 -- the maximum amount of cycles default: 200
    mutation_rate   = (5)                   -- Rate of mutatation.          default: 5
    evaluation_time = 10                    -- The eval time for each gen.  default: 10 seconds

    -- the minimum and maximum ranges for horizontal movement
    min_step        = (0.002*1000)          --l: 0.002  default: 0.001
    max_step        = (0.010*1000)          --l: 0.01   default: 0.015

    -- the minimum and maximum ranges for vertical movement
    min_vstep       = (0.001*1000)          --l: 0.001  default: 0.002
    max_vstep       = (0.020*1000)          --l: 0.02   default: 0.03

    -- phase (?)
    min_cm0         = (2*1000)
    max_cm0         = (9*1000)

    --create list of variable ranges (min/max)
    var_ranges = {}
    var_ranges[1] = {}
    var_ranges[2] = {}
    var_ranges[3] = {}
    var_ranges[4] = {}

    var_ranges[1][1] = min_step
    var_ranges[1][2] = max_step
    var_ranges[2][1] = min_vstep
    var_ranges[2][2] = max_vstep
    var_ranges[3][1] = min_cm0
    var_ranges[3][2] = max_cm0
    var_ranges[4][1] = 0 --fitness
    var_ranges[4][2] = 0 --fitness

    fitness_idx = 4 --the index where the fitness is stored for a genotype record
end

--results


-- /Global variables

--------------------------------
-- The most interesting function is "evolve()"
-- ignore most of the other functions below.
--------------------------------

-- the evolution software
function evolve()
    -- initialisation
    if (simGetScriptExecutionCount() == 0) then
        -- generate the initial pool with random genes
        population = {}
        learning_curve = {}
        population = initialize_pool(population)

        -- init the hexapod
        init_hexapod()

        -- set some values to the initial conditions
        eval_time = 0
        execution_count = 0

		local file = io.open("Hexa_LOG.csv", "w")
		file:write("Hexapod; chld. fit.; Avg fit.; Max fit:;=max(D2:D201) \n")
		file:close()
    end

    --do the evaluation if there are non-graded genotypes in the pool
    idx = check_for_ungraded_genotypes(population)
    if not (idx == (#population+1) )  then
        --evaluate if the timer is running, or start if 0
        if not (eval_time >= evaluation_time) then
            if(eval_time == 0) then --load the variables once
                update_variables_hexapod(population[idx])
            end

            --increment the execution_count
            eval_time = eval_time + simGetSimulationTimeStep()

            --let the hexapod run
            eval_hexapod()
            return
        else
            --get the fitness and store it
            population[idx][fitness_idx] = get_fitness()

            --increment the amount of times evaluations have taken place
            execution_count = execution_count + 1
            reset_hexapod()
            eval_time = 0
            print('executions;',execution_count)

			local file = io.open("Hexa_LOG.csv", "a")
			file:write(idx, ";", population[idx][fitness_idx], ";", get_average_fitness(population), ";", get_max_fitness(population), "\n")
			file:close()
            return
        end
    end

    --evolve the population for a fixed amount of times
    if not (execution_count == max_evaluations) then
        print('---')
        table.insert(learning_curve, get_average_fitness(population))

        -- select parents and let them reproduce
        child = crossover(select_parent(population), select_parent(population))

        --mutate the child with a chance
        child = mutate(child)

        --add the child to the population
        population = reduce_and_add_to_population(population, child)
        print('---')
    else
        --TODO: print the results to the console or a log file

        print(' ')
        print('printing best parent')
        print('================')
        best = print_best_parent(population)
        -- reset the fitness, the best child will be re-evaluated
        population[best][fitness_idx] = 0
        -- let this child run for a while
        evaluation_time = 99999
        print(' ')
        print('printing learning curve')
        print('================')
        --print learning curve
        for i=1, (#learning_curve) do
            print('iter:', i, 'fitn',learning_curve[i])
        end
        --end print learning curve
    end
end



-- Code of the original hexapod, most of it are as it was supplied

-- load the genotype to the hexapod, call this before evaluation
function update_variables_hexapod (genotype)
    -- variables which will evolve
    step    = genotype[1]   -- goal (max) horizontal step size
    vstep   = genotype[2]   -- goal (max) vertical step size
    cm0     = 4 --math.floor(genotype[3])
    cm1     = (4*cm0)
end

-- init function of the hexapod, also records the initial positions needed to make the reset function work
function init_hexapod()
    if (simGetScriptExecutionCount()==0) then
        -- SAVE (one time in the beginning)
        botPositionInitial = {}
        botOrientationInitial = {}
        t={simGetObjectHandle('hexapod')}
        while (#t~=0) do
           h=t[1]
           botPositionInitial[h] = simGetObjectPosition(h, -1)
           botOrientationInitial[h] = simGetObjectOrientation(h, -1)
           table.remove(t,1)
           ind=0
           child=simGetObjectChild(h,ind)
           while (child~=-1) do
              table.insert(t,child)
              ind=ind+1
              child=simGetObjectChild(h,ind)
           end
        end
        -- END SAVE
    end
    baseHandle=simGetObjectHandle('hexa_base')  -- get pointer to the base
    interModuleDelay=4      -- each leg has a delay of 4 entries in the movement table (offset in sliding window)
    xMovementTable={}       -- movement data for all directions (sliding window)
    yMovementTable={}
    zMovementTable={}
    tableLength=5*interModuleDelay+1    -- size of sliding window sufficient to keep movement data for all legs
    for i=1,tableLength,1 do            -- fill table with zeros
        table.insert(xMovementTable,0)
        table.insert(yMovementTable,0)
        table.insert(zMovementTable,0)
    end

    -- runtime variables
    phase=0 								-- phase = movement phase (finite state machine)
    r=0											-- horizontal movement position in step
    z=0											-- vertical movement position in step
    rearExtent=-0.04
    cnt=0										-- expired simulation time
    c=8
    prevSt=1

    -- KEEP CONSTANT, no evolution here
    rotation=0 									-- amount the robot rotates around its own axis (-1 to 1)
    direction=0 								-- direction the robot is travelling (0 to math.pi*2)
    forwardVel=1 								-- velocity of movement
end

-- the function called every 50ms to make the hexapod run
function eval_hexapod ()
		-- Following makes the robot first move straight, then in a circle (while keeping its orientation),
		--  then again in a circle but this time also rotating.
		-- **********************************************************************
		cnt=cnt+simGetSimulationTimeStep()

		-- Following piece of code is just to adjust for a simulationTimeStep change
		-- In that case the table lengths have to be adjusted and the interModule delays, the step sizes, etc.
		-- **********************************************************************
		st=math.floor((0.05/simGetSimulationTimeStep())+0.5)
		if (st~=prevSt) then
			c=c*(st/prevSt)
			cm0=cm0*(st/prevSt)
			cm1=cm1*(st/prevSt)
			interModuleDelay=interModuleDelay*(st/prevSt)
			otx=xMovementTable
			oty=yMovementTable
			otz=zMovementTable
			otl=tableLength
			xMovementTable={}
			yMovementTable={}
			zMovementTable={}
			tableLength=5*interModuleDelay+1
			if ((st/prevSt)>1) then
				for i=1,tableLength,1 do
					table.insert(xMovementTable,otx[math.floor((i/(st/prevSt))+0.51)])
					table.insert(yMovementTable,oty[math.floor((i/(st/prevSt))+0.51)])
					table.insert(zMovementTable,otz[math.floor((i/(st/prevSt))+0.51)])
				end
			else
				for i=1,tableLength,1 do
					table.insert(xMovementTable,otx[math.floor(i*(st/prevSt))])
					table.insert(yMovementTable,oty[math.floor(i*(st/prevSt))])
					table.insert(zMovementTable,otz[math.floor(i*(st/prevSt))])
				end
			end
			prevSt=st
		end
		-- **********************************************************************

		-- Calculate movement data at this time point
		if (phase==2) then
			r=r+step*2/st
			z=z-vstep/st
			c=c+1
			if (c>=cm0) then
				phase=0
				c=0
			end
		end
		if (phase==1) then
			r=r+step*2/st
			z=z+vstep/st
			c=c+1
			if (c>=cm0) then
				phase=2
				c=0
			end
		end
		if (phase==0) then
			r=r-step/st
			c=c+1
			if (c>=cm1) then
				phase=1
				c=0
			end
		end

		-- convert magnitude r to x-y vector according to set direction (but, we normally only go straight ahead)
		x=forwardVel*r*math.cos(direction)
		y=forwardVel*r*math.sin(direction)

		table.remove(xMovementTable,tableLength)		-- remove last item of FIFOs
		table.remove(yMovementTable,tableLength)
		table.remove(zMovementTable,tableLength)
		table.insert(xMovementTable,1,x)				-- Insert position data in FIFOs
		table.insert(yMovementTable,1,y)
		table.insert(zMovementTable,1,z)

		-- Communicate data (FIFOs) to each of the legs
		simSendData(sim_handle_tree,0,'HEXA_x',simPackFloats(xMovementTable))
		simSendData(sim_handle_tree,0,'HEXA_y',simPackFloats(yMovementTable))
		simSendData(sim_handle_tree,0,'HEXA_z',simPackFloats(zMovementTable))

		-- execute child scripts
		simHandleChildScript(sim_handle_all_except_explicit,baseHandle,interModuleDelay,rotation)
end

-- reset the hexapod to it's original start position
function reset_hexapod ()
-- RESTORE (reset to test new individual)
    -- // reset bot to its initial position
    -- // code taken from: http://www.forum.coppeliarobotics.com/viewtopic.php?f=9&t=685
    -- // apparently simSetObjectPosition() itself does not implicitly/sufficiently call simResetDynamicObject() for all the children
    t={simGetObjectHandle('hexapod')}
    while (#t~=0) do
       h=t[1]
       simResetDynamicObject(h)
       simSetObjectPosition(h, -1, botPositionInitial[h])
       simSetObjectOrientation(h, -1, botOrientationInitial[h])
       table.remove(t,1)
       ind=0
       child=simGetObjectChild(h,ind)
       while (child~=-1) do
          table.insert(t,child)
          ind=ind+1
          child=simGetObjectChild(h,ind)
       end
    end
-- END RESTORE
    --re-init makes it slightly better
    init_hexapod()
end
-- /Code of the original hexapod





-- Evolutionary software


--  generate a new genotype with random genes
--  return new_genotype
function generate_random_genotype ()
    new_genotype = {}

    new_genotype[1] = math.random(var_ranges[1][1], var_ranges[1][2])/1000  --min/max step
    new_genotype[2] = math.random(var_ranges[2][1], var_ranges[2][2])/1000  --min/max Vstep
    new_genotype[3] = math.random(var_ranges[3][1], var_ranges[3][2])/1000  --min/max cm0
    new_genotype[fitness_idx] = 0                                           -- must be 0 to tell that it isnt graded

    return new_genotype
end


--  initialize the population pool
--  return the new population filled with random genotypes
function initialize_pool(population)
    for i=1, (population_size) do
        population[i] = {}
        new_genotype = generate_random_genotype()
        population[i] = new_genotype
    end
    return population
end

--checks the pool for a ungraded genotype, returns the index if found
--returns population_size+1 if not found
function check_for_ungraded_genotypes(population)
    for i=1, (#population) do
        if population[i][fitness_idx] == 0 then
            return i
        end
    end
    return (#population+1)
end


-- select a parent using the wheel of fortune
function select_parent(population)
    total_fitness = 0.0
    rnd_point = 0.0
    fitness_sum = 0

    -- calculate the total_fitness
    for i=1, #population do
        total_fitness = total_fitness + population[i][fitness_idx] -- [i][4] -> fitness, sum fitnesses
    end

    --pick a number between 0 and the total_fitness
    rnd_point = (math.random() * total_fitness) -- select random point
    i = 1

    --get the index that corresponds to the random number
    while rnd_point > fitness_sum  do
        fitness_sum = fitness_sum + population[i][fitness_idx]
        i = i + 1
    end

    print('chosen parent id;', (i-1))
    return population[(i-1)]
end


--  crossover between two parents, the child has 50% chance to take the genotypes of either parents
--  new_child
function crossover (parent1, parent2)
    --todo misschien middelen
    new_child = {}

    for i=1, #parent1 do
        if (math.random(0,1) == 1) then
            new_child[i] = parent2[i]
        else
            new_child[i] = parent1[i]
        end
    end

    return new_child
end


--  mutate the child according to the mutation_rate.
--  mutated child
function mutate(child)
    mutated_child = {}

    for i=1, #child do
        if (math.random(0,100) < mutation_rate) then
            mutated_child[i] = math.random(var_ranges[i][1], var_ranges[i][2])/1000
        else
            mutated_child[i] = child[i]
        end
    end

    return mutated_child
end


-- return the index of the genotype with the worst fitness
function get_idx_of_worst_fitness(population)
    idx = 1

    --loop until lowest value is found
    lowest_value = 99999
    for i=1, (#population) do
        if population[i][fitness_idx] < lowest_value then
            idx = i
            lowest_value = population[i][fitness_idx]
        end
    end

    return idx
end


-- overwrite the bad gene with a new child
function reduce_and_add_to_population(population, child)
    --get a bad genotype
    idx = get_idx_of_worst_fitness(population)

    --overwrite bad genotype with the new child
    for i=1, (#child) do
        population[idx][i] = child[i]
    end

    --reset the fitness of the newly added child
    population[idx][fitness_idx] = 0
    print('new child id: ',idx)
    return population
end


--print the best parent
--also returns idx found
function print_best_parent(population)
    idx = 1

    --loop until best value is fnd
    highest_value = 0
    for i=1, (#population) do
        if population[i][fitness_idx] > highest_value then
            idx = i
            highest_value = population[i][fitness_idx]
        end
    end

    print('Best parent id:', idx, 'Fitness:', population[idx][fitness_idx])
    --print('step',population[idx][1], 'vstep',population[idx][2], 'cm0', population[idx][3])
    print('step',population[idx][1], 'vstep', population[idx][2])
    return idx
end


--
function get_average_fitness(population)
    total = 0
    for i=1, (#population) do
        total = total + population[i][fitness_idx]
    end

    --calc the average score
    average = (total/#population)

    return average
end

function get_max_fitness(population)
	idx = 1

    --loop until best value is fnd
    highest_value = -1
    for i=1, (#population) do
        if population[i][fitness_idx] > highest_value then
            idx = i
            highest_value = population[i][fitness_idx]
        end
    end

	return highest_value
end

--  get the fitness value from the hexapod, this is the distance traveled (x)
--  return fitness value
function get_fitness ()
    handle=simGetObjectHandle('hexa_body')	-- get pointer to the base
    pos = simGetObjectPosition(handle, -1)

    x = pos[1]
    -- prevent from returning 0, see also function check_for_ungraded_genotypes()
    if (x == 0) then
        return 0.001
    end

    return x
end


evolve()
-- /Evolutionary software
